NON-DEGENERATE EULERIAN FINITE ELEMENT METHOD FOR 
SOLVING PDES ON SURFACES 
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Abstract. The paper studies a method for solving elliptic partial differential equations posed on 
hypcrsurfaces in M. N , N = 2, 3. The method builds upon the formulation introduced in Bertalmio et 
al., J. Comput. Phys., 174 (2001), 759—780., where a surface equation is extended to a neighborhood 
of the surface. The resulting degenerate PDE is then solved in one dimension higher, but can be 
solved on a mesh that is unaligned to the surface. We introduce another extended formulation, which 
leads to uniformly elliptic (non-degenerate) equations in a bulk domain containing the surface. We 
apply a finite element method to solve this extended PDE and prove the convergence of finite element 
solutions restricted to the surface to the solution of the original surface problem. Several numerical 
examples illustrate the properties of the method. 
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1. Introduction. Partial differential equations posed on surfaces arise in math- 
ematical models for many natural phenomena: diffusion along grain boundaries |24j . 
lipid interactions in biomembranes [16] . and transport of surfactants on multiphase 
flow interfaces [20], as well as in many engineering and bioscience applications: vec- 
tor field visualization [13) . textures synthesis [30], brain warping [2S], fluids in lungs 
[21j among others. Thus, recently there has been a significant increase of interest in 
developing and analyzing numerical methods for PDEs on surfaces. 

The development of numerical methods based on surface triangulation started 
with the paper of Dziuk [14] . In this class of methods, the surface is approximated by 
a family of consistent regular triangulations. It is typically assumed that all vertices in 
the triangulations lie on the surface. In [15] the method from [14] was combined with 
Lagrangian surface tracking and was generalized to equations on evolving surfaces. To 
avoid surface triangulation and remeshing, another approach was taken in [7J: It was 
proposed to extend a partial differential equation from the surface to a set of positive 
Lebesgue measure in M 3 . The resulting PDE is then solved in one dimension higher, 
but can be solved on a mesh that is unaligned to the surface. If the surface evolves, 
the approach allows to avoid a Lagrangian description of the surface evolution and is 
commonly referred to as Eulerian approach [52"] . 

Despite clear advantages, the method from [7J has a number of drawbacks, see 
[101 ITS] for the careful account of pros and cons of the approach. In particular, 
the resulting bulk elliptic or parabolic PDEs are degenerate, since no diffusion acts 
in the direction normal to the surface. Setting boundary conditions in a numerical 
method for such problem is another issue. An attempt to overcome the degeneracy and 
related issues was made in [18] . where a modification of the method was introduced 
for parabolic problems. 

The method for surface equations in the present paper benefits from the modi- 
fication introduced by Greer in [T5J of the formulation from [7J. We develop a new 
extended formulation, which leads to a uniformly elliptic equations in a bulk domain 
containing the surface. The formulation preserves all advantages of the one from [TS] . 
but adds diffusion in the normal direction in a more consistent way and avoids intro- 
ducing additional parameters. Further, we consider a Galerkin (finite element) method 
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for solving the extended equation. Taking the advantage of the non-degeneracy of the 
extended formulation, we prove error estimates in the L 2 and L°° surface norms. To 
the best of our knowledge, such estimates were previously unknown for an Eulerian 
surface finite element method based on extension of PDE from the surface. 

Another Eulerian finite element method for elliptic equations posed on surfaces 
was introduced in j25j [26] . That method does not use an extension of the surface 
partial differential equation. It is instead based on a restriction (trace) of the outer 
finite element spaces to a surface. It is not the intension of this paper to compare this 
different approaches. 

The remainder of the paper is organized as follows. Section [2] collects some nec- 
essary definitions and preliminary results. In section [3J we recall the extended PDE 
approach from 7\ and the modified formulation from [18] . Further, we introduce the 
new formulation and show its well-posedness. In section]!] we consider a finite element 
method and prove error estimates. Section [5] presents the result of several numerical 
experiments that demonstrate the performance of the finite element method. Finally, 
section [6] collects some closing remarks. 

2. Preliminaries. We assume that f2 is an open subset in Mr , N = 2,3 and 
r is a connected C 2 compact hypersurface contained in fl. For a sufficiently smooth 
function g : ft — > R the tangential gradient (along T) is defined by 

Vrff = V 5 -V 3 -n r n r . (2.1) 

By Ar we denote the Laplace-Beltrami operator on T, Ar = Vr • Vr- 

This paper deals with elliptic equations posed on T. As a basic elliptic equation, 
we consider the Laplace-Beltrami problem: 

-A r u + aii = / on T, (2.2) 

with some a G L°°(r). The corresponding weak form of (2.2) reads: For given 
/ G L 2 (T) determine u G H 1 ^) such that 

J V r iiVr« + a™ds = J fvds for all v G H 1 (T). (2.3) 



For the well-possedness of (2.3), it is sufficient to assume a to be strictly positive on 



a subset of T with positive surface measure: 

A := meas s {x G T : a(x) > a } > 0, (2.4) 

with some ao > 0. In this case, the following Friedrich's type inequality [21] (see, also 
Lemma 3.1 in |25J: 

IMII^r) <C F (\\V r v\\ 2 LHr) + \\V^v\\h {r) ) VveH 1 ^) (2.5) 
holds with a constant Cf dependent of «o & n d A. 



The solution u to (2.3) is unique and satisfies u G H 2 (T), with ||u||jj2(r) < 
c ll/IU 2 (r) and a constant c independent of /, cf. [14]. We remark that the case a = 
is also covered by the analysis of the paper. In this case, the Friedrich's inequality 



(2.5) holds for all v G -ff 1 (r) with zero mean. Hence, if a = 0, we assume J r f ds = 
and look for the unique solution to (2.3) satisfying J r uds = 0. 
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Denote by ild a domain consisting of all points within a distance from T less than 
some d > 0: 

n d = {xeR N : dist(x,T) < d}. (2.6) 

Let </> : fid — > K be the signed distance function, \<fi{x)\ := dist(x,r) for all x £ fi^. 
The surface T is the zero level set of <j>: 

r = {x e R N : 0(x) = 0}. (2.7) 

We may assume <f> < on the interior of T and <f> > on the exterior. We define 
n(x) := V0(x) for all x G il^. Thus, n is the outward normal vector on T, nr = V</> 
on T, and |n(x)| = 1 for all x £ Sl d - The Hessian of (j> is denoted by H: 

H(x) = D 2 4>{x) £ M 3x3 for all x G fi d . (2.8) 

The eigenvalues of H(x) are denoted by Ki(x), «2(x), and 0. For x £ T, the eigenvalues 
Ki(x), j = 1,2, are the principal curvatures. 
We will need the orthogonal projection 

P(x) = I — n(x) <8> n(x) for all x G Qa- 

Note that the tangential gradient can be written as Vrg(x) = PVa(x) for x G T. We 
introduce a locally orthogonal coordinate system by using the projection p : Sl^ — > T: 

p(x) = x — 0(x)n(x) for all x € Sid- 

Assume that the decomposition x = p(x) + </>(x)n(x) is unique for all x G f^. We 
shall use an extension operator defined as follows. For a function u on T we define 

u e (x) := v(p(x)) for all x G fi d . (2.9) 

Thus, u e is the extension of u along normals on T, it satisfies n • Vu e = in Qd, i-e., 
u e is constant along normals to T. 

3. Extended surface PDEs. In this section, we define an extension of the 



surface PDE (2.2 1 to a neighborhood of T. Recalling (2.7 1 and Vru = PVu on T, we 



write the weak formulation of (2.2) on the zero level set of 



Vrit • Vru + auv — fv ds = / PVu • PVu + c uv — fv ds = 0. 
r J{0=o} 



The idea of [7] is to extend (2.3 1, with the help of globally defined quantities n and 



P, to every level set of 4> intersecting Sid- Find u £ Hp such that 
0=1 [ PVu • PVu + a e uv — f e v ds dr 

J-d J{4>=r} 



(3-1) 

/ (PWu-PWv + a e uv- f e v)\V(t)\dx for all v £ Hp, 
Jn d 



where 



Hp = {v £ L {Q d ) : PVu G (L 2 (Oi)) iV }- 
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The above weak formulation was shown to be well-posed in [9] . 
( 2.3 ) is embedded in (3.1| and the solution on every level set of ( 



The surface equation 
does not depend on a 



data in a neighborhood of this level set (indeed, one can consider ( 3.1 1 as a collection of 
of mutually independent surface problems on every level set of <f). Hence, restricted 
to r, smooth solution to (3.1) solves the original Laplace-Beltrami problem (2.2). 



With no ambiguity, we shall denote by u both the solutions to surface and extended 
problems. 



The corresponding strong formulation of (3.1) is 
- IV0I- 1 div |V0|PVu + a e u = 



f in n d 



(3.2) 



We note that (3.1) and (3.2) are the valid extensions of (2.3) and (2.2) if <f> is an 



arbitrary smooth level set function with V0 ^ 0, not necessarily a signed distance 
function, and a e ,f e are not necessarily constant along normal directions. If the 
boundary of the volume domain is not a level set of <j>, then (3.2) should be 



complemented with boundary conditions. This can be natural boundary conditions 



(PV«) • ngQ = on 90(j 



(3.3) 



where ngu is the outward normal vector to dQd- 

The major numerical advantage of the extended formulation is that one may apply 
standard discretization methods to solve 



( |3.2[ ) ( 3.3 1 in the volume domain Qd (e.g., 
a finite difference method on Cartesian grids) and further take the trace of computed 
solutions on T (or on a approximation of T). Numerical experiments from [7l l9l fTHl 132] 
suggest that these traces of numerical solutions are reasonably good approximations 
to the solution of the surface problem (2.2). The analysis of the method is still 



limited: Error estimates for finite element methods for (3.1) are shown in [HI ITU] . 
Error estimate in [5] is established only in the integral volume norm 



li 2 (O d ) 



IPVi 



lL 2 (n d )' 



rather than in a surface norm for T. In [lOj the first order convergence was proved in 
the surface H 1 norm, if the band width d in (2.6) is of the order of mesh size and if 



a quasi-uniform triangulation of ft is assumed. For linear elements this estimate is of 
the optimal order. 

Although numerically convenient, the extended formulation has a number of dis- 
advantages, as noted already in [7] and reviewed in [TO], [18] . The volume formulation 
(3.2) is defined in a domain in one dimension higher than the surface equation. This 
leads to involving extra degrees of freedom in numerical method. If fid is a narrow 



band around T, then handling boundary conditions (3.3) may effect the quality of the 
discrete solution. This can be an issue for grids not aligned with a level set of cf> on 
dftd- Numerical stability calls for the extension of data satisfying (2.9|; and in time- 



stepping schemes for parabolic problems, one needs the intermittent re-initialization 
of u by re-extending it from V according to (2.9). Another issue of the extended 



formulation (3.2) is that the second order term is degenerate, since no diffusion acts 
in the direction normal to level sets of <fi. Numerical solution of degenerate elliptic 
and parabolic equations is not a very well understood subject. 

An effort to overcome the degeneracy and some related issues of the approach 
from [7] was done by Greer in [TS] , where the heat equation 



du 



A r u = 0, u\ t= o = u 
l 



(3.4) 



on a stationary surface T was studied. In the method from [18], one extends (3.4 1 to 
a neighborhood of T ensuring the following properties hold: 
1. 4> is the singed distance function; 



2. Mo is extended to the neighborhood of T according to (2.9 1, i.e., constant 
alone normals; 

3. The projection P is changed to the (non-orthogonal) scaled projection 

P := (I - 0H) _1 P (3.5) 

on tangential planes of the level sets of <f>. The bulk domain Q d is assumed 
such that the modified projection P is well defined. For a smooth T, this can 
be always enured by choosing small enough d > 0. 
With the above assumptions, the solution to the extended heat equation 

3u — — 

— - (PV) • PV« = 0, u|t=o = «g in n d (3.6) 

is proved to be constant in normal directions: 

(n ■ V)u = in Q d (3.7) 

for all t > 0. 

The property (3.7) is crucial, since it allows to add diffusion in the normal di- 



rection without altering solution. Doing this, one obtains a non-degenerated elliptic 



operator. Thus, instead of (3.6) it was suggested in [18] to consider the parabolic 
problem 

du ~ ~ 

— — (PV) • PVu — c 2 n div(n ® n)Vu = 0, u| t=0 = u in Q d , (3.8) 

with a coefficient c„. For the planar case, fl d e M 2 , it was recommended to set 
c n — (1 — <Mo)i K o — K (p( x ))> K is t ne curvature of T (T is a curve in this case). For 
the case of surfaces embedded in IR 3 , there was no clear recommendation on c n . 

The above approach formally solves the problem of the degeneracy and suggests 



that equation (3.7) on d£l d is appropriate and numerically sound boundary condition. 
However, one has to define parameter c„. Moreover, the new extended formulation 
involves the Hessian H. If T is given only by an approximation, for example, as the 
zero set of a discrete level set function <f>h, then computing (an approximation to) H 
is a delicate issue, sensitive to numerical implementation. 

Below we introduce a formulation of the extended surface problem, which 'auto- 
matically' generates diffusion in the normal direction, leading to a uniformly elliptic or 
parabolic problem in il d . The finite element method and error analysis are considered 
in the section [4] The problem of the approximate evaluation of Hessian is addressed 
numerically in section [5] 

3.1. Another extension of surface PDE. For the sake of analysis, consider 



the Laplace-Beltrami equation ( 2.2 ) rather than the surface heat equation. We assume 



from now that all extensions of data from T satisfy (2.9). Consider the Laplace- 
Beltrami equation extended from r to fl d : 

-(PV)-PVu + a e u=f in fi d . (3.9) 
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To ensure that P is well-defined and equations (3.9) are well-possed, it is sufficient 



for the matrix (I — </>H) to be uniformly positive definite in f^- Therefore, assume 
f2<2 is such that 



|0(x)| = dist(x,r) < -HHMir 1 vx e n d . 



(3.10) 



One can always satisfy the above restriction by choosing the band width d small 
enough. To be more precise, from (2.5) in |11) we have the following formula for the 
eigenvalues of H: 



«i(p(x)) 



1 + ^(x)/Si(p(x)) 



for x 6 rid- 



Thus, assumption (3.10) is true if the parameter d in (2.6) satisfies 

d< (4max(|«i(x)| + |k 2 (x)|)) . 

Since T € C 2 and T is compact, the principle curvatures of T are uniformly bounded 
and d can be chosen sufficiently small positive. 



The weak formulation of the problem (3.9) reads: Find u € Hp satisfying 



PVu • PVt; + a e uv efac 



f e vdx VvGHp. 



(3.11) 



If (3.10) holds, the existence of the unique solution to (3.11) follows from the Lax- 



Milgram lemma. If the solution to (3.11) is smooth, it solves the surface problem 



(2.2) (P = P on r). Moreover, the smooth solution to (3.11) satisfies (3.7). To see 
this, apply (n • V) to equation (3.9) and use the following commutation property (see 
lemma 1 in [18] ) : 

(n • V)((PV) ■ PV) = ((PV) • PV)(n ■ V). 
Recalling (n ■ V)a e = (n ■ V)/ e = 0, we get for v n := (n • V)w 
-(PV)-PVv„ + a e v„ = in fl d . 



The uniqueness result yields (3.7) 



Note that the identity HP = PH implies 

(I - (f>B.)- 1 P = P(I - 0H)" 1 . 



Using (3.12) and P 2 = P = P J , we rewrite (3.11 1 as 



/ (I-(j)li)- 1 P\7u-(I-(j)U)- 1 \7v + a e uvdx= [ f vdx V v G H 1 
Jn d Jn d 



(3.12) 



(fid). (3.13) 



Due to relations |(I — P)Vu| 2 = |(n-Vu)n| 2 = (n-Vit) 2 , we can rewrite equality 



(3.7) for the solution to (3.13) in the form 



(I-P)Vu = 0. 
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Thanks to (3.12), it holds 



(I - ^H^V'u = P(I - 0H)- x Vii + (I - P)(I - 0H) _1 Vu 
= (I - 0H) -1 PVu + (I - ^H)- 1 ^ - P)Vu 
= (I - ^H) _1 PVu for u solving pl3||. 



We infer that the problem (3.13) can be written as follows: Find u £ i? 1 (f2 ( j) 

I {I-4>H)- 2 Vu-Vv + a e uvdx= j f e vdx for all v G H 1 (fi d ) . (3.14) 

Now we find the strong form of (3.14). To handle boundary terms arising from 
integration by part, we note that n = n^n, since the boundary of the volume domain 
fid is a level set of 4>. Furthermore, Hn = implies (I — i/)H) _1 n = n, and so 

((I - ^H) _1 V«) • n = (Vw) • ((I - 0H) _1 n) = (n • V)u. 



Thus, one can write (3.14) in the strong form: 

- div(I - cj)H)- 2 Vu + a e u = f in Q a 



(n- V)t 



on <9fL 



(3.15) 



The formulation (3.15) has the following advantages over (3.2), (3.8) and (3.9): 
The equations (3.15) are non-degenerate and uniformly elliptic, the extended problem 
has no parameters to be defined, the boundary conditions are given and consistent 
with the solution property (3.7). 



Regarding the well-posedness of (3.15) we prove the following result. 
Theorem 3.1. Assume ( |3.10| ), then it holds: 

(i) The problem (3.15) has the unique weak solution u 6 H 1 ^^, which satisfies 

\\ u \\H 1 (Q, d ) — ^ H/ e IU 2 (f2 d )i with a constant C dependent only on a, T and d; 

(ii) Additionally assume T £ C 3 , then u £ H 2 (fld) and 

IM|ff2(Q d ) < C ||/ 6 ||L2 (0£i ), 

where the constant C depends only on a, T and d. 
Proof. First we check that the bilinear form 

a(u, v) := / (I - </>H)~ 2 Vu • Vw + a e uvdx 

is continuous and coercive on H 1 (fld)- 

The assumption (3.10) yields for the spectrum of the symmetric matrices: 

sp(I - 0H) € [§, |] sp ((I - 0H)- 2 ) G [f,4] for any x £ fl d . 

Therefore, it holds 

5l|V«||i» (n< ) < / (I-^H)- 2 V U -V U dx, (3.16) 

f (I-0H)- 2 V U -V«dx<4||V U || L 2 (Od) ||V U || L 2 (f2d) . (3.17) 



Estimates ( |3.17 1 and ||a e ||L°°(o ci ) = IMlL°°(r) imply the continuity estimate 



Define 



a(u,v) < 4||VM|| L 2 (0ti) ||Vw|| L 2 (s:2d) + \\a e \\ L ^ (nd) \\u\\ L 2 {nd) \\v\\ L 2 {nd) 
< (4+\\a\\ L ~ { r))\\u\\ H i{n d )\\v\\m(n d )- 



M(x) := (1 -d(x)Ki(x))(l -d(x)K 2 (x)), xe fl d . 



From (2.20), (2.23) in jTT] we have /i(x)dx = drds(p(x)), for x £ 17^ where dx is the 
measure in fid, ds the surface measure on T, and r the local coordinate at x £ T in 



the normal direction. Using (3.101, we get \ < /x(x) < | for all x 6 fi^. From this 



and relations (2.9 1 and (2.4), we infer that a e is strictly positive on a subset of fid 



with positive measure: 



A := meas x {x £ fid : a e (x) > a } > „d meas s {x e T : a(x) > a } > 0, 

9 



with «q > 0. Hence, similar to the surface case in ( 2.5 ), the Friedrich's type inequality 



\\v\\l HQd) < C F (\\Vv\\ 2 L2m + \\V^vf L2(Qa) ) V u £ H 1 ^) (3.18) 
holds with a constant Cf dependent of ao and A. 



Inequalities (3.16) and (3.18) imply the ellipticity of the bilinear form: a{u,u) > 
c\\u\\ 2 for all u £ i? 1 (Sl^) , where the constant c depends only on Cf from (3.18). 
Therefore, part (i) of the theorem follows from the Lax-Milgram lemma. 

To check part (ii) of the theorem, we note that T £ C 3 yields 
and dfld £ C 3 . Therefore, the entries of the 'diffusion' matrix (I — 0H)~ 2 are in C 1 
and a £ L°°(T) => a e £ L°°(f2d). This smoothness of the data is sufficient for the 
elliptic problem to be i? 2 -regular [2j [22] and the result follows. 
□ 

Remark 1 . Theorem |3.1| shows one theoretical advantage of the new extended 



£ C 3 , see Ql], 
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formulation (3.15) over (3.9) and (3.2): If the data is smooth, then the Agmon- 
Douglis-Nirenberg regularity theory immediately applies. In particular, for r £ C 3 , 
we have u £ H 2 {Q d ) and the trace theorem, see, e.g., [22], yields u\r £ i? 1 (r). This 



enables one to consider the trace of u as the weak solution to (2.2) 



4. Finite element method. Let V £ C 2 and fix a domain such that the 



band width d satisfies (3.10). Assume T is a consistent division (triangulation) of 
fid into tetrahedra elements. We call a triangulation of fid exact if lJ Tg7 -T = fid- 
Since dfld coincides with isolines of the distance function <f», the boundary of fid is 
curvilinear: dfld £ C 2 . Hence, exact triangulations of fid may be constructed only in 
certain cases using isogeometric elements [3] or mapped (blending) finite elements [33] . 
In a general case, we define the domain: 

fifc := IJ T < 

TET 

which approximates fid- 

Furthermore, in some applications the surface T may not be known explicitly, but 
given only approximately as, for example, the zero level set of a finite element distance 
function ^h- In this case, instead of the Hessian H = V 2 4> one has to use a discrete 



Hessian H/j H, which is obtained from ^ by any of discrete Hessian recovery 
methods, see e.g. [5J|3T]. We assume that <j>h and satisfy condition (3.10). 

Let Vh C if 1 (fi^,) be a space of finite element functions. The finite element 
method reads: Find Uh G 14 satisfying 

/ ( (I - <j> h H h )- 2 Vu h ) ■ Vv h + a e u h u fc dx = / /Sftdx Vu^ G V h . (4.1) 

We analyse the method (4.1 ) below in the special case of O/j = f2<£, = 0, and 
Hft = H. Numerical experiments in the next section test the method when non of 
these assumptions hold. 

Since the diffusion tensor (I — </)H)~ 2 is uniform positive definite and bounded, 
we immediately obtain the following optimal convergence result, e.g., [8]: 

Theorem 4.1. Let f2ft, = f^ 7 fa = (j), and H/> = H. Assume u and solve 
problems (3.15) and (4.1), respectively. Then it holds 

\u-v h \\ H i( Qi ), 



\u - u h \\ H i {nd) < C inf 



where the constant C may depend only on a, T and d. 

Theorem |4 . 1 1 and the trace theorem yield the simple error estimate on the surface: 



\\u-u h \\ L -2 (T) <C inf \\u - v h \\ H i (nd) . 



(4.2) 



However, the above estimate of the error in surface L 2 norm is not optimal and can 
be improved. The improved estimate is given by Theorem |4.5| To show it, we need 
several preparatory results. 
Denote 



h = sup diam(T) 

TGT 



and assume that Vh is such that 



inf ||u - V h \\ H irQ d \ < C a h\\v\\ H 2r Ud ), Vv£H 2 (fL). 

Vh&Vh 



(4.3) 



The L 2 -convergence estimate for the finite element method for the extended problem 
is given in the next theorem. 

Theorem 4.2. Let fi/,. = fid, fa — <t>> H/j = H, and T G C 3 . Assume u and Uh 
solve problems (3.15) and (4.1), respectively. Then it holds 

Uh\\L,i(ci d ) < Ch inf \\u-v h \\ H un d ), 

Vh£V h 



where the constant C may depend only on a, T, d, and constant C a from (4.3). 

Proof. The assumption r G C 3 ensures that the differential problem (3.15) is 
H 2 -regular. Since flh = ^d, <t>h = <t>, = H, the discrete problem (4.1) is the plain 
Galerkin method. Hence, the result follows from the standard duality argument, see, 
e.g., 0. □ 

The result below is found for example in Theorem 1.4.3.1 of [19] (note that the 
unit simplex has uniformly Lipschitz boundary). 

Lemma 4.3. Let T be the unit simplex (triangle or tetrahedra) in R , N = 2,3. 
Then there exists an extension operator E : iJ x (T) — > H 1 (M. N ) such that 



|£v|liji(K«) < C \\v\\ Hl{f) V v G H 1 ^). 



(4.4) 



For a tetrahedron (triangle) T denote by p{T) the diameter of the inscribed ball. 
Let 7r be the set of all tetrahedra intersected by T. Denote 



(3 = sup diam(T)/p(T). 

TeTr 



(4.5) 



We assume that tetrahedra (triangles) in 7r are shape- regular, i.e., (3 is not too big. 
We need the following technical lemma. 

Lemma 4.4. Let T G 7r- Denote h = diam(T) and K = T n F, then it holds 



\M L 2 ( K) <C( h 2 \Mlht) +h*\\Vv\\ L 2 {T )), 



(4.6) 



where the constant C may depend only on F and the constant (3 from (4.5|. 

Proof. Fhe proof adopts the 'flattening' argument from 12J, § 3.4. Fhe proof 
below is given for the three-dimensional case: T is a surface in M. 3 . All arguments 
remain valid with obvious modifications, if T in a curve in M 2 . We may assume 
that the curvilinear element K has non-zero 2D measure. Let T be the reference 
unit tetrahedron in R 3 , let ip : T — > T be an affine mapping with ||V<p|| < ch and 
||(V(p) _1 || < ch -1 . Here and in the rest of the proof, c denotes a generic constant, 
which may depend only on /? and T, but does not depend on T. We next recall 
that because r is a C 2 surface, there exists a C 2 chart $ with uniformly bounded 
derivatives, and for which <& _1 has uniformly bounded derivatives, which maps an 
0(l)-neighborhood N of K in M 3 to M 3 and which has the property that T Pi N lies 
in a plane. It is not difficult to extend to all of M 3 so that the resulting extension 
has bounded derivatives, has a bounded inverse, and flattens an 0(l)-neighborhood 
of K. We then define a corresponding flattening map for the reference space by 
<f> = ip^ 1 o $ o ip. It is easy to check that then $ and are also uniformly bounded 
in C 2 , and ^(^^-(K)) is flat. Denote by P a plane in M 3 containing the flattened 
surface element <i>((p _1 (A')). 

We need the following trace inequality ([T], Theorem 7.58): 



'\l 2 (p) < 



£3) 



(4.7) 



Define v on T by v = v o ip. Given T € T, recalling the definition of the extension 



operator E from Lemma 4.3 and trace inequality (4.7 1, we then compute 



'L 2 (K\ 



<C 

=c 
<c 
<c 
<c 
<c 
<c 



\L2( v -i(K)) 



Ev\ 



L2(<p-i(K)) 



Evo $ 



lL3 ( <E, (y -i (K))) 



Evo $ _1 || i 2( ] 
Evo 
Ev\ 



(4.8) 



ff 1 



\h 1 (t)- 



Applying a scaling argument yields 



\H X (T) 



<c(^ 3/2 ||«l| z3(T) +/i- 1 /*||V?;|| £a( D). 
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(4.9) 



Estimates (4.8) and (4.9) prove the lemma. □ 

Summing up the estimate ( |4.6[ ) over all elements from 7r, we get for v € H 1 (Qd) 



\v\\h { r)<C(hJ n ]T \\v\\l HT) + /w H Vw l 

TeTr TGTr 

<C(h m l\\v\\ 2 L2(nd) +h max \\Vv\\ 2 L2(nd) ), 



i 2 (T)) 



(4.10) 



with /i m ; n ( max ) = min(max)T e 7- r diam(T). For the next theorem, let us assume h max < 

C ^min • 

Theorem 4.5. Lei Qh — f^, 4>h — 4>- H/j = H, and Y £ C 3 . Assume u and Uh 



solve problems (2.2) and (4.1), respectively. Then it holds 



\u - u h \\ L 2 (r) < Ch? inf \\u-v h \\ H i/ a ^, 



where the constant C may depend only on a, Y, d, constants C a from (4.3) , and f3 
from (4.5 I. 



Proof. The result of the theorem follows from Theorems 4.1 4.2 and (4.10). □ 



The error estimate from Theorem 4.5 is an improvement of (4.2), but still half- 



order suboptimal. The difficulty in proving the optimal estimate is that the analysis 
of the extended finite element problem in 1 (Sl^) and L 2 (VLd) norms gives little in- 
formation of normal derivatives of the error, i.e. how accurate Uh satisfies (3.7). 



Approaching optimal estimates of the error on Y is possible applying the well- 
known results on interior maximum-norm estimates for finite element methods for el- 
liptic problems |27j . This, however, requires some further restrictions on mesh and Vh- 
To be precise, assume a quasi-uniform triangulation of Qh and let dist(<97r, dtth) > 
C\d > c h, with ci independent of h and large enough constant Cq. Let Vh be the 
space of Pk finite element functions, k > 1. Now we apply Theorem 1.1 from [2"7] 
(in terms of |27j we take the 'basic' domain Qq = 7r and the intermediate domain 
ild = Qh)- We obtain 



\u-Uh\\L°°(V) < \\u-u h \\ L oa (Tr) 

<C[{lndh' 1 ) r min \\u - Vh\\ L °°(n h ) + d~% \\u - u h \\ L 2 {Qh) ] 

\ v h £V h J 



(4.11) 



for u and Uh solving (3.15) and (4.1), respectively. Here r = 1 for k — 1 and r = 

for k > 2. 

Now we want to combine the result in (4.11) with L 2 volume estimate from 
Theorem |4.2| To do this, we have to assume tth 



fid, which means curvilinear 
elements on the boundary of O/j. Although certain types of curvilinear elements are 
allowed by the analysis of [27], we avoid further assumptions on elements touching 
boundary, but simply separate from the boundary: Consider fl' h = {T £ T : TO 
dttd — 0}- Assume £l' h consists only of shape-regular tetrahedra (triangles). The 
restriction of finite element functions from Vh on VL' h is denoted by Vh(£l' h ). We 
assume Vh{Q' h ) is the space of P k elements. If d is fixed and h is sufficiently small, 
then dist(<97r, d£l' h ) > c\d > cqH, with c\ independent of h and large enough constant 
cq. Hence, the result in (4.11) holds with Qh replaced by fl' h and we can combine it 



with the estimate from Theorem |4.2| Thus, we proved the following theorem. 

Theorem 4.6. Let£l h = il d , 4>h = <j>, = H, T e C 3 , d is fixed such that (13.101 



holds, h is sufficiently small, the triangulation of Q' h is quasi-uniform and consists of 
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Fig. 5.1. The triangulation of after one step of refinement for the 2D example and the 
visualization of the computed solution of the extended problem. 



tetrahedra (triangles), and V/,.(0^) is the space of P k elements, k > 1. A 



ssume u ana 



Uh solve problems (2.2) and (4.1), respectively. Then it holds 



\u - u h \\L°°(r) 



< C (ln/i 



mm 

Vh£V h 



U-V h \\ L oo ( Q h) + 



h inf 



Vh\\m(n d ) 



r = 1 for k = 1 and r = for k > 2. The constant C may depend only on a, T, d, 



constants C a from (4.3), and (3 from (4.5) 



As an example, assume u is sufficiently smooth and Vh is piecewise linear finite 
element space (mapped piecewise linear near dfld). Then Theorem |4 . 6| yields optimal 
order convergence result (up to logarithmic term): \\u — u/i||l°°(T) = 0(h 2 In h^ 1 ). 

5. Numerical examples. In this section, we present results of several numerical 
experiments. We start with the example of the Laplace-Beltrami problem (2.2) on 
a unit circle in M 2 with a known solution so that we are able to calculate the error 
between the continuous and discrete solutions. We set a — 1 and consider 

u(r, (f>) = cos(5</>) 

in polar coordinates, similar to the Example 5.1 from [10] . 

For d = 0.05, we build conforming quasi-uniform triangulation of Qd and apply 
the regular refinement process. The grid is always aligned with the boundary of fid 
so that dflh is an 0(h 2 ) approximation of dfld- The grid after one step of refinement 
in the upper right part of dilh is shown in Figure [5~T| (left). 



The solution computed on grid level 3 is shown in Figure 5.1 (right). The vi- 



sualized solution looks constant in normal direction, as expected for solution of the 
extended problem. Next, we compute the finite element error restricted to the sur- 
face. For evaluating this error, we consider piecewise linear approximations of F and 
evaluate the errors alon g t hese piecewise linear s urfac es. To assess the estimates 
given by Theorems 4.5 and 4.6 we show in Table 5.1 the surface L 2 and C-norms 
of the errors. We clearly see the second order of convergence in both norms, when 
the Hessian of the distance function is taken in (4.1) exactly. We also experiment 
with approximate choice of 4>h and H/j in (4.1 ). In this case, 4>h is a piecewise linear 
Lagrange interpolant to <f>, and is a piecewise linear continuous tensor-function 
recovered from (f)^ by the variation method [4J. Compared to the exact choice, the 
finite element errors are somewhat larger, although the convergence rates stay close 
to the second order. The discrete problems were solved using the BCG method with 
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Table 5.1 

Norms of the errors for the example of sphere with exact and approximate Hessian. n,j is fixed 
with d = 0.1. # Iter is the number of preconditioned BCG iterations. 





h 


#d.o.f. 


L^-norm 


Order 


C-norm 


Order 


# Iter. 




0.0417 


610 


0.318E-02 




0.345E-02 




13 




0.0208 


2058 


0.662E-03 


2.26 


0.148E-02 


1.22 


28 


H 


0.0104 


7351 


0.179E-03 


1.89 


0.308E-03 


2.26 


60 




0.0052 


27954 


0.409E-04 


2.13 


0.812E-04 


1.92 


142 




0.0026 


109576 


0.983E-05 


2.06 


0.195E-04 


2.06 


325 




0.0417 


610 


0.449E-02 




0.451E-02 




13 




0.0208 


2058 


0.147E-02 


1.61 


0.182E-02 


1.31 


28 




0.0104 


7351 


0.390E-03 


1.91 


0.420E-03 


2.16 


63 




0.0052 


27954 


0.124E-03 


1.65 


0.160E-03 


1.39 


137 




0.0026 


109576 


0.325E-04 


1.93 


0.425E-04 


1.91 


297 



Table 5.2 

Norms of the errors for the example of sphere with exact and approximate Hessian, fl^ is fixed 
with d = 0.1. 





#d.o.f. 


L^-norm 


Order 


C-norm 


Order 


# Iter. 




1026 


0.6085E-01 




0.9033E-01 




9 


H 


8547 


0.1503E-01 


1.98 


0.1523E-01 


2.52 


23 




63632 


0.3990E-02 


1.98 


0.3971E-02 


2.01 


47 




1026 


0.8095E-01 




0.1032E+00 




9 


H h 


8547 


0.2144E-01 


1.88 


0.1909E-01 


2.39 


23 




63632 


0.5114E-02 


2.14 


0.4529E-02 


2.15 


46 



the ILU2 preconditioner [53] to a relative tolerance of 10 -9 . The iterations numbers 
are shown in the right column of the table. 

As the next test problem, we consider the Laplacc-Bcltrami equation on the unit 
sphere: 

— Arit + u — / on F, 

with T = {x € R 3 | 1 1 x 1 1 2 = 1}- The source term / is taken such that the solution is 
given by 

12 

u(x) = pT|3 (3xjx 2 - x 2 ) , x = (xi, x 2 , x 3 ) e fi. 

Note that u and / are constant along normals at T. 

For different values of the domain width parameter d, we build conforming sub- 
divisions of fid into regular-shaped tetrahedra using the software package ANI3D [3J . 
The grid is aligned with the boundary of fid so that dflh is an 0(h 2 ) approximation 
of dfld- The resulting discrete problems are again solved by the BCG method with 
the ILU2 preconditioner to a relative tolerance of 10~ 9 . 

In Tables |5.2| |5.3| we show the norms of surface errors for the computed finite 
element solutions and the number of preconditioned BCG iterations. The surface 
errors were computed using the piecewise planar approximations of T by T^, where 
T/i is the zero level set of the piecewise linear Lagrange interpolant to the distance 
function of T. Table [5T2| shows the error norms for the case of a fixed domain and 
a sequence of discretizations. The formal convergence order p was computed as 



p = 3 log (erri/err 2 ) / log ((#d.o./.)i/(#d.o./.) 2 ) 
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Fig. 5.2. The visualization of solution on the sphere and the cutaway of the volume grid in 
for d = {0.1, 0.2, 0.4}. 



Table 5.3 

Dependence of error norms on the domain width d for the example with sphere. 



d 


#d.o.f. 


L^-norm 


C-norm 


# Iter. 


0.4 


80442 


0.7700E-02 


0.6986E-02 


46 


0.2 


34305 


0.9389E-02 


0.9560E-02 


12 


0.1 


13560 


0.9579E-02 


0.1025E-01 


35 



We clearly observe the second order convergence both when the exact distance func- 



tion (j) and the Hessian H was used in (4.1) and when (j> and H were replaced by 



discrete (fih and H^. As in the two-dimension case, <ph is a piecewise linear Lagrange 
interpolant to <fi an d is a piecewise linear continuous tensor-function recovered 
from 4>h by the variation method. In Table |5.3[ the results are shown for the case 
when the domain fl^ shrinks towards the surface, while the mesh size h was approxi- 
mately the same for all three meshes. 



The Figure 5.2 visualizes the solutions for various widths of the volume domains 
fid- We see that the discrete solutions tend to be constant in the normal direction to 
the surface. 

We repeat the previous experiment, but now with a torus instead of the unit 
sphere. Let F = {x e Q | r 2 = x\ + (y/x\ + x 2 ~ R ) 2 }- We take R = 1 and r = °- 6 - 
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Fig. 5.3. The visualization of solution on the torus and the cutaway of the volume grid in 
for d = 0.1. 

Table 5.4 

Norms of the errors for the example of torus with exact and approximate Hessian. f2(j is fixed 
with d = 0.1. 





#d.o.f. 


L 2 -norm 


Order 


C-norm 


Order 


# Iter. 




26257 


0.7826E-01 




0.1405E+00 




18 


H 


174021 


0.2843E-01 


1.61 


0.8400E-01 


0.82 


42 




1511742 


0.7780E-02 


1.80 


0.1077E-01 


2.85 


98 




26257 


0.1144E+00 




0.1708E+00 




20 




174021 


0.7680E-01 


0.63 


0.1569E+00 


0.13 


13 




1511742 


0.6888E-01 


0.15 


0.9679E-01 


0.67 


94 



In the coordinate system (p, <f>, 9) , with 




cos 4> cos < 
p I sin cos ( 
sinf? 



the p-direction is normal to T, |£ _L T for x e T. The following solution u and 
corresponding right-hand side / are constant in the normal direction: 

u(x) = sin(30) cos(36* + </>), 
/(x) = r- 2 (9 sin(3</>) cos(3<9 + </>)) 

-(R + r cos(9)- 2 (-10 sin(30) cos(36< + <j>) - 6 cos(3<£) sin(36» + 0)) 

- (r(R + r cos^))- 1 (3 sin(0) sin(30) sin(36» + 0)). 



(5.1) 



The surface norms of approximation errors for the example of torus are given in 
Table |5.4| The solution is visualized in Figure |5.3| Again, when the exact Hessian is 
used, the convergence order is close to the second one. However, when the exact Hes- 
sian is replaced by the recovered Hessian, the convergence significantly deteriorates. 
This is opposite to the example with the sphere. A closer inspection reveals that the 
error is concentrated in the proximity of the inner ring of the torus, where the Gauss 



curvature is negative (see Figure 5.4 left). Next, we look on the err or o f the discrete 

(EL =1 (H 



Hessian recovery: 
that the error |H 



|H- 



H,, 



The Figure 



5.4 



right, shows 
of the torus. 



is large at the same region, near the inner ring 
At this part of fid the Hessian is indefinite, it has non-zero values of different signs 
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--0.108 



- -0.267 



- -0.214 



--0.161 



—0.0547 



—0.0016 



- 0.051 



- 0.104 



- 0.157 




- 0.01 



Fig. 5.4. Left: The finite element method error, with approximate Hessian. Right: The error 
between discrete and approximate Hessian. 

6. Conclusions. We studied a formulation and a finite element method for el- 
liptic partial differential equation posed on hypersurfaces in R N , N = 2,3. The 
formulation uses an extension of the equation off the surface to a volume domain 
containing the surface. Unlike the original method from [7], the extension introduced 
in the paper results in uniformly elliptic problems in the volume domain. This en- 
ables a straightforward application of standard discretization techniques and put the 
problem into a well-established framework for analysis of elliptic PDEs, including 
numerical analysis. For the standard Galerkin finite element method we proved new 
convergence estimates in the surface L? and L°° norm. Optimal convergence of the PI 
finite clement method was demonstrated numerically. The method, however, requires 
a reasonably good approximation of the Hessian for the signed distance function to 
the surface. 
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